  module met_preproc_mod
    use spark_cubed_sphere
    use params_mod
    use dycore_test
    use static_mod
    use state_mod
    use ghost_interp_mod
    use topo_mod
    use vertical_coord_mod
    use phy_operator_mod
    use fiona
  
    implicit none
    
    type met_field_type
      integer :: nlon
      integer :: nlat
      integer :: nlev
      integer :: nlev_soil
      
      real(r_kind) :: min_lon, max_lon
      real(r_kind) :: min_lat, max_lat
      real(r_kind) :: dlon, dlat
      
      real(r_kind), dimension(:  ), allocatable :: lon      
      real(r_kind), dimension(:  ), allocatable :: lat      
      real(r_kind), dimension(:,:), allocatable :: lon_2d   
      real(r_kind), dimension(:,:), allocatable :: lat_2d   
      real(r_kind), dimension(:  ), allocatable :: mlev     
      real(r_kind), dimension(:,:), allocatable :: mlev_2d  
      real(r_kind), dimension(:  ), allocatable :: plev     
      real(r_kind), dimension(:,:), allocatable :: plev_2d  
      real(r_kind), dimension(:  ), allocatable :: zlev     
      real(r_kind), dimension(:,:), allocatable :: zlev_2d  
                                                            
      real(r_kind), dimension(:,:,:), allocatable :: gh     
      real(r_kind), dimension(:,:,:), allocatable :: rho    
      real(r_kind), dimension(:,:,:), allocatable :: p      
      real(r_kind), dimension(:,:,:), allocatable :: t      
      real(r_kind), dimension(:,:,:), allocatable :: pt
      real(r_kind), dimension(:,:,:), allocatable :: u      
      real(r_kind), dimension(:,:,:), allocatable :: v      
      real(r_kind), dimension(:,:,:), allocatable :: w      
      real(r_kind), dimension(:,:,:), allocatable :: sh     
      real(r_kind), dimension(:,:,:), allocatable :: rh     
      real(r_kind), dimension(:,:,:), allocatable :: mr     
                                                            
      real(r_kind), dimension(:,:  ), allocatable :: ps     
      real(r_kind), dimension(:,:  ), allocatable :: u10m   
      real(r_kind), dimension(:,:  ), allocatable :: v10m   
      real(r_kind), dimension(:,:  ), allocatable :: t2m    
      real(r_kind), dimension(:,:  ), allocatable :: d2m    
      real(r_kind), dimension(:,:  ), allocatable :: q2m     
      real(r_kind), dimension(:,:  ), allocatable :: mr2m   
      real(r_kind), dimension(:,:,:), allocatable :: slc    
      real(r_kind), dimension(:,:,:), allocatable :: smc    
      
      real(r_kind), dimension(:,:  ), allocatable :: zs
      
      contains
      procedure :: clear => met_field_clear
      final :: met_final
    end type
    
    type interim_field_type
      ! 1: npoints
      ! 2: nx
      ! 3: ny
      ! 4: nlev
      !                       1 2 3 4
      real(r_kind), dimension(:,:,:,:), allocatable :: gh     
      real(r_kind), dimension(:,:,:,:), allocatable :: rhod   
      real(r_kind), dimension(:,:,:,:), allocatable :: rhom   
      real(r_kind), dimension(:,:,:,:), allocatable :: rho    
      real(r_kind), dimension(:,:,:,:), allocatable :: p      
      real(r_kind), dimension(:,:,:,:), allocatable :: t      
      real(r_kind), dimension(:,:,:,:), allocatable :: pt
      real(r_kind), dimension(:,:,:,:), allocatable :: u      
      real(r_kind), dimension(:,:,:,:), allocatable :: v      
      real(r_kind), dimension(:,:,:,:), allocatable :: w      
      real(r_kind), dimension(:,:,:,:), allocatable :: sh     
      real(r_kind), dimension(:,:,:,:), allocatable :: rh     
      real(r_kind), dimension(:,:,:,:), allocatable :: mr     
                                                            
      real(r_kind), dimension(:,:  ,:), allocatable :: ps     
      real(r_kind), dimension(:,:  ,:), allocatable :: u10m   
      real(r_kind), dimension(:,:  ,:), allocatable :: v10m   
      real(r_kind), dimension(:,:  ,:), allocatable :: t2m    
      real(r_kind), dimension(:,:  ,:), allocatable :: d2m    
      real(r_kind), dimension(:,:  ,:), allocatable :: q2m     
      real(r_kind), dimension(:,:  ,:), allocatable :: mr2m   
      real(r_kind), dimension(:,:,:,:), allocatable :: slc    
      real(r_kind), dimension(:,:,:,:), allocatable :: smc    
      
      real(r_kind), dimension(:,:  ,:), allocatable :: zs
      
      contains
      procedure :: clear => interim_field_clear
      final :: interim_final
    end type
    
    type(met_field_type    )              :: met
    type(interim_field_type), allocatable :: interim(:)
    type(interim_field_type), allocatable :: fields(:)
    
  contains
    
    subroutine met_get_data_era_interim(high_field_file, surface_field_file, invariant_field_file, itime)
      character(*), intent(in) :: high_field_file
      character(*), intent(in) :: surface_field_file
      character(*), intent(in) :: invariant_field_file
      integer     , intent(in) :: itime
      
      real(r_kind), dimension(:,:,:  ), allocatable :: tmp3d
      real(r_kind), dimension(:,:,:,:), allocatable :: tmp4d
      
      real(r_kind) :: R, Cp, Cv, alpha, mc ! moisture correction coef
      
      integer i,j,k
      
      call log_notice('Get initial data from ERA-interim file.', pid=proc%id)
      
      call fiona_init()
      call fiona_open_dataset('high_field', file_path=trim(high_field_file))
      call fiona_get_dim('high_field', 'longitude', size=met%nlon)
      call fiona_get_dim('high_field', 'latitude' , size=met%nlat)
      call fiona_get_dim('high_field', 'level'    , size=met%nlev)
      
      allocate(met%lon  (met%nlon))
      allocate(met%lat  (met%nlat))
      allocate(met%plev (met%nlev))
      allocate(met%gh   (met%nlon, met%nlat, met%nlev))
      allocate(met%t    (met%nlon, met%nlat, met%nlev))
      allocate(met%u    (met%nlon, met%nlat, met%nlev))
      allocate(met%v    (met%nlon, met%nlat, met%nlev))
      allocate(met%w    (met%nlon, met%nlat, met%nlev))
      allocate(met%sh   (met%nlon, met%nlat, met%nlev))
      allocate(met%mr   (met%nlon, met%nlat, met%nlev))
      allocate(met%pt   (met%nlon, met%nlat, met%nlev))
      allocate(met%p    (met%nlon, met%nlat, met%nlev))
      
      allocate( tmp3d(met%nlon, met%nlat,           1) )
      allocate( tmp4d(met%nlon, met%nlat, met%nlev, 1) )
      
      call fiona_start_input('high_field')
      call fiona_input('high_field', 'longitude', met%lon )
      call fiona_input('high_field', 'latitude' , met%lat )
      call fiona_input('high_field', 'level'    , met%plev)
      call fiona_input('high_field', 'z'        , tmp4d, start=(/1,1,1,itime/), count=(/met%nlon,met%nlat,met%nlev,1/))
      met%gh = tmp4d(:,:,:,1)
      call fiona_input('high_field', 't'        , tmp4d, start=(/1,1,1,itime/), count=(/met%nlon,met%nlat,met%nlev,1/))
      met%t = tmp4d(:,:,:,1)
      call fiona_input('high_field', 'u'        , tmp4d, start=(/1,1,1,itime/), count=(/met%nlon,met%nlat,met%nlev,1/))
      met%u = tmp4d(:,:,:,1)
      call fiona_input('high_field', 'v'        , tmp4d, start=(/1,1,1,itime/), count=(/met%nlon,met%nlat,met%nlev,1/))
      met%v = tmp4d(:,:,:,1)
      call fiona_input('high_field', 'w'        , tmp4d, start=(/1,1,1,itime/), count=(/met%nlon,met%nlat,met%nlev,1/))
      met%w = tmp4d(:,:,:,1)
      call fiona_input('high_field', 'q'        , tmp4d, start=(/1,1,1,itime/), count=(/met%nlon,met%nlat,met%nlev,1/))
      met%sh = tmp4d(:,:,:,1)
      call fiona_end_input('high_field')
      
      call fiona_open_dataset('surface_field', file_path=trim(surface_field_file))
      allocate(met%t2m (met%nlon, met%nlat))
      allocate(met%d2m (met%nlon, met%nlat))
      allocate(met%u10m(met%nlon, met%nlat))
      allocate(met%v10m(met%nlon, met%nlat))
      allocate(met%ps  (met%nlon, met%nlat))
      allocate(met%mr2m(met%nlon, met%nlat))
      
      call fiona_start_input('surface_field')
      call fiona_input('surface_field', 't2m', tmp3d, start=(/1,1,itime/), count=(/met%nlon,met%nlat,1/))
      met%t2m = tmp3d(:,:,1)
      call fiona_input('surface_field', 'd2m', tmp3d, start=(/1,1,itime/), count=(/met%nlon,met%nlat,1/))
      met%d2m = tmp3d(:,:,1)
      call fiona_input('surface_field', 'u10', tmp3d, start=(/1,1,itime/), count=(/met%nlon,met%nlat,1/))
      met%u10m = tmp3d(:,:,1)
      call fiona_input('surface_field', 'v10', tmp3d, start=(/1,1,itime/), count=(/met%nlon,met%nlat,1/))
      met%v10m = tmp3d(:,:,1)
      call fiona_input('surface_field', 'sp' , tmp3d, start=(/1,1,itime/), count=(/met%nlon,met%nlat,1/))
      met%ps = tmp3d(:,:,1)
      call fiona_end_input('surface_field')
      
      call fiona_open_dataset('invariant_field', file_path=trim(invariant_field_file))
      allocate(met%zs(met%nlon, met%nlat))
      
      call fiona_start_input('invariant_field')
      call fiona_input('invariant_field', 'z', met%zs)
      call fiona_end_input('invariant_field')
      
      met%lon = met%lon * D2R
      met%lat = met%lat * D2R
      
      met%max_lon = maxval(met%lon)
      met%max_lat = maxval(met%lat)
      
      met%min_lon = minval(met%lon)
      met%min_lat = minval(met%lat)
      
      met%dlon = met%lon(2) - met%lon(1)
      met%dlat = met%lat(2) - met%lat(1)
      
      if(met%dlat<0)then
        call flipud(met%lat)
        met%dlat = abs(met%dlat)
        do k = 1, met%nlev
          do i = 1, met%nlon
            call flipud( met%gh(i,:,k) )
            call flipud( met%t (i,:,k) )
            call flipud( met%u (i,:,k) )
            call flipud( met%v (i,:,k) )
            call flipud( met%w (i,:,k) )
            call flipud( met%sh(i,:,k) )
            call flipud( met%mr(i,:,k) )
            call flipud( met%p (i,:,k) )
          enddo
        enddo
      endif
      
      if( met%plev(2) - met%plev(1) > 0 )then
        call flipud( met%plev )
        do j = 1, met%nlat
          do i = 1, met%nlon
            call flipud( met%gh(i,j,:) )
            call flipud( met%t (i,j,:) )
            call flipud( met%u (i,j,:) )
            call flipud( met%v (i,j,:) )
            call flipud( met%w (i,j,:) )
            call flipud( met%sh(i,j,:) )
            call flipud( met%mr(i,j,:) )
            call flipud( met%p (i,j,:) )
          enddo
        enddo
      endif
      
      ! Convert data
      met%plev = met%plev * 100
      met%mr = met%sh / ( 1._r_kind - met%sh )
      !met%mr = 0

      do k = 1, met%nlev
        do j = 1, met%nlat
          do i = 1, met%nlon
            met%p(i,j,k) = met%plev(k)
            call calc_thermal_parameters(met%mr(i,j,k), R, Cp, Cv, met%sh(i,j,k), alpha=alpha, mc=mc)
            met%pt(i,j,k) = met%t(i,j,k) * ( p0 / met%p(i,j,k) )**( R / Cp )
          enddo
        enddo
      enddo
      
      do j = 1, met%nlat
        do i = 1,met%nlon
          met%mr2m(i,j) = calc_mr_from_dew_point_temperature(met%d2m(i,j), met%t2m(i,j), met%ps(i,j))
        enddo
      enddo
      
      print*,'min/max value of lon  ', minval(met%lon  ), maxval(met%lon  )
      print*,'min/max value of lat  ', minval(met%lat  ), maxval(met%lat  )
      print*,'min/max value of plev ', minval(met%plev ), maxval(met%plev )
      print*,'min/max value of gh   ', minval(met%gh   ), maxval(met%gh   )
      print*,'min/max value of t    ', minval(met%t    ), maxval(met%t    )
      print*,'min/max value of u    ', minval(met%u    ), maxval(met%u    )
      print*,'min/max value of v    ', minval(met%v    ), maxval(met%v    )
      print*,'min/max value of w    ', minval(met%w    ), maxval(met%w    )
      print*,'min/max value of sh   ', minval(met%sh   ), maxval(met%sh   )
      print*,'min/max value of mr   ', minval(met%mr   ), maxval(met%mr   )
      print*,'min/max value of pt   ', minval(met%pt   ), maxval(met%pt   )
      print*,'min/max value of mr2m ', minval(met%mr2m ), maxval(met%mr2m )
      print*,'min/max value of u10m ', minval(met%u10m ), maxval(met%u10m )
      print*,'min/max value of v10m ', minval(met%v10m ), maxval(met%v10m )
      print*,'min/max value of ps   ', minval(met%ps   ), maxval(met%ps   )
      print*,'min/max value of zs   ', minval(met%zs   ), maxval(met%zs   )
      print*,''
      
    end subroutine met_get_data_era_interim
    
    subroutine met_interp_fields(state, static)
      type(state_type ), intent(inout) :: state
      type(static_type), intent(in   ) :: static
      
      integer idom, ims, ime, jms, jme, kms, kme, ids, ide, jds, jde, kds, kde
      integer ndom, npt, pc, pvs, pve, pes(6), pee(6), pqs, pqe, pgs, pge
      integer i, j, k, p
      integer is, ie, js, je, ks, ke
      real(r_kind) :: lonL, lonR, latB, latT
      real(r_kind) :: R, Cp, Cv, alpha, mc ! moisture correction coef
      real(r_kind) :: dp, rhog, dh
      real(r_kind) :: iAh(2,2)
      
      call log_notice('Interpolating meterological fields to interim grid.', pid=proc%id)
      
      associate (mesh => state%mesh)
        call mesh%get_params(npt=npt, pc=pc, pvs=pvs, pve=pve, pes=pes, pee=pee, pqs=pqs, pqe=pqe, pgs=pgs, pge=pge)
        allocate( interim(mesh%ndom) )
        allocate( fields (mesh%ndom) )
        
        do idom = 1, mesh%ndom
          if(mesh%domain(idom)%initialized)then
            call mesh%get_params(idom=idom,                                            &
                                 ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
                                 ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde)
            
            allocate( interim(idom)%p    (pqs:pqe,ims:ime,jms:jme,0:met%nlev) )
            allocate( interim(idom)%gh   (pqs:pqe,ims:ime,jms:jme,0:met%nlev) )
            allocate( interim(idom)%t    (pqs:pqe,ims:ime,jms:jme,0:met%nlev) )
            allocate( interim(idom)%u    (pqs:pqe,ims:ime,jms:jme,0:met%nlev) )
            allocate( interim(idom)%v    (pqs:pqe,ims:ime,jms:jme,0:met%nlev) )
            allocate( interim(idom)%w    (pqs:pqe,ims:ime,jms:jme,0:met%nlev) )
            allocate( interim(idom)%mr   (pqs:pqe,ims:ime,jms:jme,0:met%nlev) )
            allocate( interim(idom)%pt   (pqs:pqe,ims:ime,jms:jme,0:met%nlev) )
            allocate( interim(idom)%ps   (pqs:pqe,ims:ime,jms:jme           ) )
            allocate( interim(idom)%u10m (pqs:pqe,ims:ime,jms:jme           ) )
            allocate( interim(idom)%v10m (pqs:pqe,ims:ime,jms:jme           ) )
            allocate( interim(idom)%t2m  (pqs:pqe,ims:ime,jms:jme           ) )
            allocate( interim(idom)%mr2m (pqs:pqe,ims:ime,jms:jme           ) )
            allocate( interim(idom)%zs   (pqs:pqe,ims:ime,jms:jme           ) )
            
            associate( lon   => mesh%domain(idom)%lon,&
                       lat   => mesh%domain(idom)%lat,&
                       press => interim(idom)%p,      &
                       gh    => interim(idom)%gh,     &
                       t     => interim(idom)%t,      &
                       u     => interim(idom)%u,      &
                       v     => interim(idom)%v,      &
                       w     => interim(idom)%w,      &
                       mr    => interim(idom)%mr,     &
                       pt    => interim(idom)%pt,     &
                       ps    => interim(idom)%ps,     &
                       u10m  => interim(idom)%u10m,   &
                       v10m  => interim(idom)%v10m,   &
                       t2m   => interim(idom)%t2m,    &
                       mr2m  => interim(idom)%mr2m,   &
                       zs    => interim(idom)%zs      )
              ! interpolate high fields
              !$OMP PARALLEL
              !$OMP DO PRIVATE(is, ie, js, je, lonL, lonR, latB, latT) COLLAPSE(4)
              do k = 1, met%nlev
                do j = jms, jme
                  do i = ims, ime
                    do p = pqs, pqe
                      !call ll2ij(is, ie, js, je, lon(p,i,j), lat(p,i,j), met%min_lon, met%max_lon, met%min_lat, met%max_lat, met%dlon, met%dlat, met%nlon, met%nlat)
                      call ll2ij_directly(is, ie, js, je, lon(p,i,j), lat(p,i,j), met%lon, met%lat)
                      lonL = met%lon(is)
                      lonR = merge( met%lon(ie) + 2 * pi, met%lon(ie), ie < is)
                      latB = met%lat(js)
                      latT = met%lat(je)
                      gh(p,i,j,k) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%gh(is,js,k), met%gh(ie,js,k), met%gh(is,je,k), met%gh(ie,je,k))
                      t (p,i,j,k) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%t (is,js,k), met%t (ie,js,k), met%t (is,je,k), met%t (ie,je,k))
                      u (p,i,j,k) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%u (is,js,k), met%u (ie,js,k), met%u (is,je,k), met%u (ie,je,k))
                      v (p,i,j,k) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%v (is,js,k), met%v (ie,js,k), met%v (is,je,k), met%v (ie,je,k))
                      w (p,i,j,k) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%w (is,js,k), met%w (ie,js,k), met%w (is,je,k), met%w (ie,je,k))
                      mr(p,i,j,k) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%mr(is,js,k), met%mr(ie,js,k), met%mr(is,je,k), met%mr(ie,je,k))
                      pt(p,i,j,k) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%pt(is,js,k), met%pt(ie,js,k), met%pt(is,je,k), met%pt(ie,je,k))
                      
                      press(p,i,j,k) = met%plev(k)
                    enddo
                  enddo
                enddo
              enddo
              !$OMP END DO NOWAIT
              
              ! interpolate surface fields
              !$OMP DO PRIVATE(is, ie, js, je, lonL, lonR, latB, latT) COLLAPSE(3)
              do j = jms, jme
                do i = ims, ime
                  do p = pqs, pqe
                    !call ll2ij(is, ie, js, je, lon(p,i,j), lat(p,i,j), met%min_lon, met%max_lon, met%min_lat, met%max_lat, met%dlon, met%dlat, met%nlon, met%nlat)
                    call ll2ij_directly(is, ie, js, je, lon(p,i,j), lat(p,i,j), met%lon, met%lat)
                    lonL = met%lon(is)
                    lonR = merge( met%lon(ie) + 2 * pi, met%lon(ie), ie < is)
                    latB = met%lat(js)
                    latT = met%lat(je)
                    ps  (p,i,j) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%ps  (is,js), met%ps  (ie,js), met%ps  (is,je), met%ps  (ie,je))
                    u10m(p,i,j) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%u10m(is,js), met%u10m(ie,js), met%u10m(is,je), met%u10m(ie,je))
                    v10m(p,i,j) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%v10m(is,js), met%v10m(ie,js), met%v10m(is,je), met%v10m(ie,je))
                    t2m (p,i,j) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%t2m (is,js), met%t2m (ie,js), met%t2m (is,je), met%t2m (ie,je))
                    mr2m(p,i,j) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%mr2m(is,js), met%mr2m(ie,js), met%mr2m(is,je), met%mr2m(ie,je))
                    zs  (p,i,j) = bilinear_interp(lon(p,i,j), lat(p,i,j), lonL, lonR, latB, latT, met%zs  (is,js), met%zs  (ie,js), met%zs  (is,je), met%zs  (ie,je))
                  enddo
                enddo
              enddo
              !$OMP END DO
              !$OMP END PARALLEL
              
              ! Insert surface data at bottom
              do j = jms, jme
                do i = ims, ime
                  do p = pqs, pqe
                    press(p,i,j,0) = 0
                    !if( ps(p,i,j) > press(p,i,j,1) )then
                    !  dp   = press(p,i,j,1) - ps(p,i,j)
                    !  rhog = press(p,i,j,1) / ( Rd * t2m (p,i,j) * ( 1._r_kind + eqv * mr2m(p,i,j) / ( 1._r_kind + mr2m(p,i,j) ) ) ) * earth_gravity
                    !  dh   = - dp / rhog
                    !  gh   (p,i,j,0) = min( static%domain(idom)%zsQ(p,i,j) * earth_gravity, gh(p,i,j,1) - dh, zs(p,i,j) )
                    !  press(p,i,j,0) = ps  (p,i,j)
                    !  t    (p,i,j,0) = t2m (p,i,j)
                    !  u    (p,i,j,0) = u10m(p,i,j)
                    !  v    (p,i,j,0) = v10m(p,i,j)
                    !  w    (p,i,j,0) = 0
                    !  mr   (p,i,j,0) = mr2m(p,i,j)
                    !endif
                  enddo
                enddo
              enddo
              
              !print*,'min/max value of gh  ', minval(gh  ), maxval(gh  )
              !print*,'min/max value of t   ', minval(t   ), maxval(t   )
              !print*,'min/max value of u   ', minval(u   ), maxval(u   )
              !print*,'min/max value of v   ', minval(v   ), maxval(v   )
              !print*,'min/max value of w   ', minval(w   ), maxval(w   )
              !print*,'min/max value of mr  ', minval(mr  ), maxval(mr  )
              !print*,'min/max value of ps  ', minval(ps  ), maxval(ps  )
              !print*,'min/max value of u10m', minval(u10m), maxval(u10m)
              !print*,'min/max value of v10m', minval(v10m), maxval(v10m)
              !print*,'min/max value of t2m ', minval(t2m ), maxval(t2m )
              !print*,'min/max value of mr2m', minval(mr2m), maxval(mr2m)
              !print*,'min/max value of zs  ', minval(zs  ), maxval(zs  )
              !print*,''
            end associate
              
            call log_notice('Interpolating meterological fields to computational grid.', pid=proc%id)
      
            ! Vertical interpolation
            allocate( fields(idom)%gh   (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%pt   (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%t    (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%p    (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%u    (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%v    (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%w    (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%sh   (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%mr   (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%rho  (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%rhod (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            allocate( fields(idom)%rhom (pqs:pqe,ims:ime,jms:jme,kms:kme) )
            
            associate( lon    => mesh%domain(idom)%lon,&
                       lat    => mesh%domain(idom)%lat,&
                       z      => mesh%domain(idom)%z,  &
                       Jab    => mesh%domain(idom)%Jab,&
                       p_src  => interim(idom)%p,      &
                       gh_src => interim(idom)%gh,     &
                       pt_src => interim(idom)%pt,     &
                       t_src  => interim(idom)%t,      &
                       u_src  => interim(idom)%u,      &
                       v_src  => interim(idom)%v,      &
                       w_src  => interim(idom)%w,      &
                       mr_src => interim(idom)%mr,     & 
                       gh     => fields(idom)%gh,      &
                       pt     => fields(idom)%pt,      &
                       t      => fields(idom)%t,       &
                       press  => fields(idom)%p,       &
                       u      => fields(idom)%u,       &
                       v      => fields(idom)%v,       &
                       w      => fields(idom)%w,       &
                       sh     => fields(idom)%sh,      &
                       mr     => fields(idom)%mr,      &
                       rho    => fields(idom)%rho,     &
                       rhod   => fields(idom)%rhod,    &
                       rhom   => fields(idom)%rhom,    &
                       q      => state%domain(idom)%q  )
              !$OMP PARALLEL DO PRIVATE(p, ks, ke, k, Cp, Cv, mc, R, alpha, iAh) COLLAPSE(2)
              do j = jms, jme
                do i = ims, ime
                  do p = pqs, pqe
                    if( p_src(p,i,j,0) <= 1.e-10 )then
                      ks = 1
                      ke = met%nlev
                    else
                      ks = 0
                      ke = met%nlev
                    endif
                    
                    ! interpolate pressure
                    call interp_1d( gh_src(p,i,j,ks:ke)/earth_gravity, log(p_src(p,i,j,ks:ke)), z(p,i,j,kds:kde), press(p,i,j,kds:kde) )
                    press(p,i,j,kds:kde) = exp(press(p,i,j,kds:kde))
                    
                    ! interpolate fields
                    call interp_1d( log(p_src(p,i,j,ks:ke)), t_src (p,i,j,ks:ke), log(press(p,i,j,kds:kde)), t (p,i,j,kds:kde) )
                    call interp_1d( log(p_src(p,i,j,ks:ke)), u_src (p,i,j,ks:ke), log(press(p,i,j,kds:kde)), u (p,i,j,kds:kde) )
                    call interp_1d( log(p_src(p,i,j,ks:ke)), v_src (p,i,j,ks:ke), log(press(p,i,j,kds:kde)), v (p,i,j,kds:kde) )
                    call interp_1d( log(p_src(p,i,j,ks:ke)), w_src (p,i,j,ks:ke), log(press(p,i,j,kds:kde)), w (p,i,j,kds:kde) )
                    call interp_1d( log(p_src(p,i,j,ks:ke)), mr_src(p,i,j,ks:ke), log(press(p,i,j,kds:kde)), mr(p,i,j,kds:kde) )
                    call interp_1d( log(p_src(p,i,j,ks:ke)), pt_src(p,i,j,ks:ke), log(press(p,i,j,kds:kde)), pt(p,i,j,kds:kde) )
                    
                    do k = kds, kde
                      call calc_thermal_parameters(mr(p,i,j,k), R, Cp, Cv, sh(p,i,j,k), alpha, mc)
                      t    (p,i,j,k) = pt(p,i,j,k) * ( press(p,i,j,k) / p0 )**( R / Cp )
                      !pt   (p,i,j,k) = t(p,i,j,k) * ( p0 / press(p,i,j,k) )**( R / Cp )
                      rhod (p,i,j,k) = press(p,i,j,k) / ( Rd * t(p,i,j,k) )
                      rhom (p,i,j,k) = rhod(p,i,j,k) * mr(p,i,j,k)
                      rho  (p,i,j,k) = rhod(p,i,j,k) + rhom(p,i,j,k)
                      
                      ! Convert wind to computational space
                      iAh = mesh%get_dom_iAh(idom,p,i,j,k)
                      call spherev_to_contrav(iAh, u(p,i,j,k), v(p,i,j,k), u(p,i,j,k), v(p,i,j,k))
                    enddo
                  enddo
                  
                  do k = kds, kde
                    q(i,j,k,1) = mesh%domain(idom)%cell_quad( Jab(pqs:pqe,i,j,k) * rhod(pqs:pqe,i,j,k)                     ) / Jab(pc,i,j,k)
                    q(i,j,k,2) = mesh%domain(idom)%cell_quad( Jab(pqs:pqe,i,j,k) * rho (pqs:pqe,i,j,k) * u (pqs:pqe,i,j,k) ) / Jab(pc,i,j,k)
                    q(i,j,k,3) = mesh%domain(idom)%cell_quad( Jab(pqs:pqe,i,j,k) * rho (pqs:pqe,i,j,k) * v (pqs:pqe,i,j,k) ) / Jab(pc,i,j,k)
                    q(i,j,k,4) = mesh%domain(idom)%cell_quad( Jab(pqs:pqe,i,j,k) * rho (pqs:pqe,i,j,k) * w (pqs:pqe,i,j,k) ) / Jab(pc,i,j,k)
                    q(i,j,k,5) = mesh%domain(idom)%cell_quad( Jab(pqs:pqe,i,j,k) * rhod(pqs:pqe,i,j,k) * pt(pqs:pqe,i,j,k) ) / Jab(pc,i,j,k)
                    q(i,j,k,6) = mesh%domain(idom)%cell_quad( Jab(pqs:pqe,i,j,k) * rhom(pqs:pqe,i,j,k)                     ) / Jab(pc,i,j,k)
                    !q(i,j,k,4) = 0
                    !q(i,j,k,6) = 0
                  enddo
                enddo
              enddo
              !$OMP END PARALLEL DO
              
              print*,'min/max value of q1', minval(q(ids:ide,jds:jde,kds:kde,1)), maxval(q(ids:ide,jds:jde,kds:kde,1))
              print*,'min/max value of q2', minval(q(ids:ide,jds:jde,kds:kde,2)), maxval(q(ids:ide,jds:jde,kds:kde,2))
              print*,'min/max value of q3', minval(q(ids:ide,jds:jde,kds:kde,3)), maxval(q(ids:ide,jds:jde,kds:kde,3))
              print*,'min/max value of q4', minval(q(ids:ide,jds:jde,kds:kde,4)), maxval(q(ids:ide,jds:jde,kds:kde,4))
              print*,'min/max value of q5', minval(q(ids:ide,jds:jde,kds:kde,5)), maxval(q(ids:ide,jds:jde,kds:kde,5))
              print*,'min/max value of q6', minval(q(ids:ide,jds:jde,kds:kde,6)), maxval(q(ids:ide,jds:jde,kds:kde,6))
            end associate
          endif
        enddo
      end associate
      
    end subroutine met_interp_fields
    
    subroutine met_field_clear(this)
      class(met_field_type), intent(inout) :: this
      
      if(allocated(met%lon    )) deallocate(met%lon    )
      if(allocated(met%lat    )) deallocate(met%lat    )
      if(allocated(met%lon_2d )) deallocate(met%lon_2d )
      if(allocated(met%lat_2d )) deallocate(met%lat_2d )
      if(allocated(met%mlev   )) deallocate(met%mlev   )
      if(allocated(met%mlev_2d)) deallocate(met%mlev_2d)
      if(allocated(met%plev   )) deallocate(met%plev   )
      if(allocated(met%plev_2d)) deallocate(met%plev_2d)
      if(allocated(met%zlev   )) deallocate(met%zlev   )
      if(allocated(met%zlev_2d)) deallocate(met%zlev_2d)
      
      if(allocated(met%gh )) deallocate(met%gh )
      if(allocated(met%rho)) deallocate(met%rho)
      if(allocated(met%p  )) deallocate(met%p  )
      if(allocated(met%t  )) deallocate(met%t  )
      if(allocated(met%u  )) deallocate(met%u  )
      if(allocated(met%v  )) deallocate(met%v  )
      if(allocated(met%w  )) deallocate(met%w  )
      if(allocated(met%sh )) deallocate(met%sh )
      if(allocated(met%rh )) deallocate(met%rh )
      if(allocated(met%mr )) deallocate(met%mr )
      if(allocated(met%pt )) deallocate(met%pt )
      
      if(allocated(met%ps  )) deallocate(met%ps  )
      if(allocated(met%u10m)) deallocate(met%u10m)
      if(allocated(met%v10m)) deallocate(met%v10m)
      if(allocated(met%t2m )) deallocate(met%t2m )
      if(allocated(met%d2m )) deallocate(met%d2m )
      if(allocated(met%q2m )) deallocate(met%q2m )
      if(allocated(met%mr2m)) deallocate(met%mr2m)
      if(allocated(met%slc )) deallocate(met%slc )
      if(allocated(met%smc )) deallocate(met%smc )
      
      if(allocated(met%zs  )) deallocate(met%zs  )
      
    end subroutine met_field_clear
    
    subroutine met_final(this)
      type(met_field_type), intent(inout) :: this

      call this%clear()
    
    end subroutine met_final
    
    subroutine interim_field_clear(this)
      class(interim_field_type), intent(inout) :: this
      
        if(allocated(this%gh  )) deallocate(this%gh  )
        if(allocated(this%rho )) deallocate(this%rho )
        if(allocated(this%p   )) deallocate(this%p   )
        if(allocated(this%t   )) deallocate(this%t   )
        if(allocated(this%u   )) deallocate(this%u   )
        if(allocated(this%v   )) deallocate(this%v   )
        if(allocated(this%w   )) deallocate(this%w   )
        if(allocated(this%sh  )) deallocate(this%sh  )
        if(allocated(this%rh  )) deallocate(this%rh  )
        if(allocated(this%mr  )) deallocate(this%mr  )
        if(allocated(this%pt  )) deallocate(this%pt  )
        if(allocated(this%ps  )) deallocate(this%ps  )
        if(allocated(this%u10m)) deallocate(this%u10m)
        if(allocated(this%v10m)) deallocate(this%v10m)
        if(allocated(this%t2m )) deallocate(this%t2m )
        if(allocated(this%d2m )) deallocate(this%d2m )
        if(allocated(this%q2m )) deallocate(this%q2m )
        if(allocated(this%mr2m)) deallocate(this%mr2m)
        if(allocated(this%slc )) deallocate(this%slc )
        if(allocated(this%smc )) deallocate(this%smc )
        if(allocated(this%zs  )) deallocate(this%zs  )
      
    end subroutine interim_field_clear
    
    subroutine interim_final(this)
      type(interim_field_type), intent(inout) :: this

      call this%clear()
    
    end subroutine interim_final
    
    subroutine ll2ij(is, ie, js, je, lon, lat, min_lon, max_lon, min_lat, max_lat, dlon, dlat, nlon, nlat)
      integer        , intent(out) :: is, ie, js, je
      real   (r_kind), intent(in ) :: lon, lat
      real   (r_kind), intent(in ) :: min_lon, max_lon
      real   (r_kind), intent(in ) :: min_lat, max_lat
      real   (r_kind), intent(in ) :: dlon, dlat
      integer(i_kind), intent(in ) :: nlon, nlat
      
      is = ( lon - min_lon ) / dlon + 1
      ie = merge( is + 1, 1, is/=nlon )
      
      js = ( lat - min_lat ) / dlat + 1
      
      if( js==nlat ) js = nlat - 1
      
      je = js + 1
      
    end subroutine ll2ij
    
    subroutine ll2ij_directly(is, ie, js, je, lon, lat, lon_src, lat_src)
      integer        , intent(out) :: is, ie, js, je
      real   (r_kind), intent(in ) :: lon, lat
      real   (r_kind), intent(in ) :: lon_src(:), lat_src(:)
      
      integer ires(1), jres(1)
      
      ires = maxloc(lon_src,lon_src<lon)
      jres = maxloc(lat_src,lat_src<lat)
      
      is = merge( ires(1), lbound(lon_src,1), ires(1)/=0 )
      js = merge( jres(1), lbound(lat_src,1), jres(1)/=0 )
      
      ires = minloc(lon_src,lon_src>lon)
      jres = minloc(lat_src,lat_src>lat)
      
      ie = merge( ires(1), ubound(lon_src,1), ires(1)/=0 )
      je = merge( jres(1), ubound(lat_src,1), jres(1)/=0 )
      
      if( je < js )then
        je = merge( js+1, ubound(lat_src,1), js < ubound(lat_src,1) )
      endif
      
    end subroutine ll2ij_directly
    
  end module met_preproc_mod
